#based on oct.genmatch.pscore
# orignally from
# lapo:/home/adiamond/midwest/DW/forJas/oct.genmatch.pscore

set.seed(2423045)
rm(list=ls())
library(Matching)

load("DW.cps1.re74.RData")
foo <- DW.cps1.re74
#to save space
rm(DW.cps1.re74)

X <-   as.matrix(cbind( foo$age,
                        foo$educ,
                        foo$black,
                        foo$hispan,
                        foo$married,
                        foo$nodegree,
                        I(foo$re74/1000),
                        #I(as.real(foo$re74 == 0)),
                        I(foo$re75/1000)))# ,
                        #I(as.real(foo$re75 == 0))))
                 
BalanceMat <- as.matrix(cbind(
                     foo$age,
                     foo$educ,
                     foo$black,
                     foo$hispan,
                     foo$married,
                     foo$nodegree,
                     I(foo$re74/1000),
                     I(foo$re75/1000),
                     I((foo$re74/1000)^2),
                     I((foo$re75/1000)^2),
                         
                     I((foo$age)^2),
                     I((foo$educ)^2),
                     I((foo$black)^2),
                     I((foo$hispan)^2),
                     I((foo$married)^2),
                     I((foo$nodegree)^2),
                         
                     I(foo$age*foo$educ),
                     I(foo$age*foo$black),
                     I(foo$age*foo$hispan),
                     I(foo$age*foo$married),    
                     I(foo$age*foo$nodegree),
                     I(foo$age*(foo$re74/1000)),
                     I(foo$age*(foo$re75/1000)),

                     I(foo$educ*foo$black),
                     I(foo$educ*foo$hispan),
                     I(foo$educ*foo$married),    
                     I(foo$educ*foo$nodegree),
                     I(foo$educ*(foo$re74/1000)),
                     I(foo$educ*(foo$re75/1000)),

                     I(foo$black*foo$hispan),
                     I(foo$black*foo$married),
                     I(foo$black*foo$nodegree),
                     I(foo$black*(foo$re74/1000)),
                     I(foo$black*(foo$re75/1000)),

                     I(foo$hispan*foo$married),
                     I(foo$hispan*foo$nodegree),
                     I(foo$hispan*(foo$re74/1000)),
                     I(foo$hispan*(foo$re75/1000)),

                     I(foo$married*foo$nodegree),
                     I(foo$married*(foo$re74/1000)),
                     I(foo$married*(foo$re75/1000)),
                                                
                     I(foo$nodegree*(foo$re74/1000)),
                     I(foo$nodegree*(foo$re75/1000)),                        

                     I((foo$re74/1000)*(foo$re75/1000))                                                 
                         ) )


### ORTHOGONALIZE TO PROPENSITY SCORE

# From Review of Econ and Statistics
model.A = foo$treat~I(foo$age) + I(foo$age^2) + I(foo$age^3) + I(foo$educ) + I(foo$educ^2) + I(foo$married) + I(foo$nodegree) + I(foo$black) + I(foo$hispan) + I(foo$re74) + I(foo$re75) + I(foo$re74 == 0) + I(foo$re75 == 0) + I(I(foo$re74)*I(foo$educ))

pscores.A <- glm(model.A, family = binomial)

X2 <- X
for (i in 1:ncol(X2))
  {
   lm1 = lm(X2[,i]~pscores.A$linear.pred)
   X2[,i] = lm1$residual
   }


# VARIABLES WE MATCH ON BEGIN WITH PROPENSITY SCORE
orthoX2.plus.pscore = cbind(pscores.A$linear.pred, X2)

orthoX2.plus.pscore[,1] <- orthoX2.plus.pscore[,1] -
                           mean(orthoX2.plus.pscore[,1])


# NORMALIZE ALL COVARS BY STANDARD DEVIATION
for (i in 1:(dim(orthoX2.plus.pscore)[2]))
  {
    orthoX2.plus.pscore[,i] <-
      orthoX2.plus.pscore[,i]/sqrt(var(orthoX2.plus.pscore[,i]))
  }

# SET BOUNDS ON GENMATCH DOMAIN (OPTIONAL)
lower <- c(95, rep(0, dim(orthoX2.plus.pscore)[2] - 1))
upper <- c(105, rep(1000, dim(orthoX2.plus.pscore)[2] - 1))

Domains <- as.matrix(cbind(lower, upper))

set.seed(12345)

source("~/R/cluster1.R")

GM.out <-
  GenMatch(Tr = foo$treat,
           X = orthoX2.plus.pscore,
           BalanceMatrix = BalanceMat,
           pop.size = 5000,
#           pop.size = 50,           
           max.generations = 25,
           wait.generations = 20,
           print.level = 3,
           cluster=cl)

#benchmark=$1794
mout <- Match(Y=foo$re78,
           Tr = foo$treat,
           X = orthoX2.plus.pscore,
           Weight.matrix=GM.out)
summary(mout)
